
rm(list=ls())

library(list)

load("data.RData")

###
###  Fig 1 - Descriptive Plots by Provinces
###

data$Q4Afac <- factor(data$Q4A, levels = c(0, 1, 2, 3, 4))
data$Q4Bfac <- factor(data$Q4B, levels = c(0, 1, 2, 3, 4))
data$Q4Cfac <- factor(data$Q4C, levels = c(0, 1, 2, 3, 4))

data$Q5.1fac <- factor(data$Q5.1, levels = c(1:5, 98, 99)) 
data$Q5.1Afac <- factor(data$Q5.1A, levels = c(1:5, 98, 99)) 
data$Q5.1Bfac <- factor(data$Q5.1B, levels = c(1:5, 98, 99)) 
data$Q5.2fac <- factor(data$Q5.2, levels = c(1:5, 98, 99)) 
data$Q5.2Afac <- factor(data$Q5.2A, levels = c(1:5, 98, 99)) 
data$Q5.2Bfac <- factor(data$Q5.2B, levels = c(1:5, 98, 99)) 
data$Q5.3fac <- factor(data$Q5.3, levels = c(1:5, 98, 99)) 
data$Q5.3Afac <- factor(data$Q5.3A, levels = c(1:5, 98, 99)) 
data$Q5.3Bfac <- factor(data$Q5.3B, levels = c(1:5, 98, 99)) 
data$Q5.4fac <- factor(data$Q5.4, levels = c(1:5, 98, 99)) 
data$Q5.4Afac <- factor(data$Q5.4A, levels = c(1:5, 98, 99)) 
data$Q5.4Bfac <- factor(data$Q5.4B, levels = c(1:5, 98, 99)) 


yaxis <- c("Control", "ISAF")

col.seq <- c(grey(seq(0.9^2.2, 0.3^2.2, length = 5)^(1/2.2)), rep(rgb(0,0,0,0), 2)) 

n.samp <- c(nrow(data[data$list.treat != "taliban", ]), rep(NA, 5))

pdf("figs/CompareEndorseListProvince.pdf", height = 5.5, width = 6.6)
par(mfrow = c(6,5), las = 1, mar=c(0,0.5,2,0.25), oma=c(4,6,2,0), cex = 0.55)

## list

ynames <- c("Control", "ISAF")

barplot(t(rbind(matrix(prop.table(table(data$Q4Afac)[1:5]), nrow = 1),
                matrix(prop.table(table(data$Q4Cfac)[1:5]), nrow = 1)
                )),
        horiz = TRUE, xaxt = "n", names.arg = yaxis, main = "", col = col.seq,
        font.main = 2, cex.main = 1.25)

mtext("List Experiment", side = 3, outer = TRUE, at = 0.1, cex = 0.65, line = 0.5, font = 2)

barplot(t(rbind(matrix(prop.table(table(data$Q5.1fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31), main = "Direct Elections"
        , font.main = 1, cex.main = 1.1) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.1fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.1Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.2fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31), main = "Prison Reform", font.main = 1, cex.main = 1.1) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.2fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.2Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.3fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31), main = "Independent Election\nCommission", font.main = 1, cex.main = 1.1) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.3fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.3Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

barplot(t(rbind(matrix(prop.table(table(data$Q5.4fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
        angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31), main = "Anti-Corruption\nReform", font.main = 1, cex.main = 1.1) 

barplot(t(rbind(matrix(prop.table(table(data$Q5.4fac)[c(1:5, 6, 7)]), nrow = 1),
                matrix(prop.table(table(data$Q5.4Afac)[c(1:5, 6, 7)]), nrow = 1))),
        horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
        col = col.seq)

par(xpd = NA)
lines(c(-4.95, 1), rep(-0.35, 2), lty = "dashed")
lines(rep(-3.3, 2), c(4.7, -27.5), lty = "dashed")

for (i in 1:length(levels(data$M5))) {

  dat <- data[data$M5 == levels(data$M5)[i], ]

  n.samp[i+1] <- nrow(dat[dat$list.treat != "taliban", ])
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q4Afac)[1:5]), nrow = 1),
                  matrix(prop.table(table(dat$Q4Cfac)[1:5]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = yaxis, main = "", col = col.seq,
          font.main = 2, cex.main = 1.25)
     
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.1fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.1fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.1Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.2fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.2fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.2Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.3fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.3fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.3Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.4fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, density = c(rep(-1, 5), rep(25, 2)),
          angle = c(rep(-45, 6), 45), lwd = 1.25, col = rgb(.31, .31, .31)) 
  
  barplot(t(rbind(matrix(prop.table(table(dat$Q5.4fac)[c(1:5, 6, 7)]), nrow = 1),
                  matrix(prop.table(table(dat$Q5.4Afac)[c(1:5, 6, 7)]), nrow = 1))),
          horiz = TRUE, xaxt = "n", names.arg = NULL, add = TRUE,
          col = col.seq)
  
}

mtext("Endorsement Experiments", side = 3, outer = TRUE, at = 0.6, cex = 0.65, line = 0.5, font = 2)

par(las = 0)
labels <- c("Overall", levels(data$M5))
for (i in 1:length(labels)) {
  mtext(labels[i], outer = TRUE, side = 2, line = 4.5, at = 0.9 - 0.165*(i-1), cex = 0.75)
  mtext(paste("(N = ", n.samp[i],")", sep=""), outer = TRUE, side = 2, line = 3.5, at = 0.9 - 0.165*(i-1), cex = 0.5)
}

par(xpd=NA)
legend(-3.7, 0, c("0 items", "1", "2", "3", "4"), 
       fill = grey(seq(0.9^2.2, 0.3^2.2, length = 5)^(1/2.2)), xjust=0.5, bty = "n", cex = 1, ncol = 2)

par(xpd=NA)
legend(-1.1, 0, c("Strongly\ndisagree","Disagree","Indifferent","Agree","Strongly\nagree","Don't know", "Refused"),
       density = c(rep(-1, 5), rep(40, 2)), angle = c(rep(-45, 6), 45),
       fill = c(col.seq[1:5], rep(rgb(.31, .31, .31), 2)), horiz=T, xjust=0.5,
       x.intersp = 0.25, bty = "n")

dev.off()

cor.dat <- cbind(Y[data$list.treat=="control" | data$list.treat == "ISAF", ],
             indivcov[data$list.treat=="control" | data$list.treat == "ISAF", ],
                   data[data$list.treat=="control" | data$list.treat == "ISAF", c("list.y","list.treat")])
cor.dat$treat <- ifelse(cor.dat$list.treat == "ISAF", 1, 0)
names(cor.dat)[1:4] <- paste("Q", 1:4, sep="")

data$treat <- data$list.treat 

data <- data[data$treat == "ISAF" | data$treat == "control", ]

data$treat.int <- as.numeric(data$treat == "ISAF")

data$dkref <- ifelse(is.na(data$Q5.1 == 98 | data$Q5.1 == 99 | data$Q5.1A == 98 | data$Q5.1A == 99 | data$Q5.1B == 98 | data$Q5.1B == 99 |
         data$Q5.2 == 98 | data$Q5.2 == 99 | data$Q5.2A == 98 | data$Q5.2A == 99 | data$Q5.2B == 98 | data$Q5.2B == 99 |
         data$Q5.3 == 98 | data$Q5.3 == 99 | data$Q5.3A == 98 | data$Q5.3A == 99 | data$Q5.3B == 98 | data$Q5.3B == 99 |
         data$Q5.4 == 98 | data$Q5.4 == 99 | data$Q5.4A == 98 | data$Q5.4A == 99 | data$Q5.4B == 98 | data$Q5.4B == 99)==F,1,0)

data.mean <- data[(data$Q5.1 <= 5 | is.na(data$Q5.1)) & (data$Q5.1A <= 5 | is.na(data$Q5.1A)) & (data$Q5.1B <= 5 | is.na(data$Q5.1B)) &
             (data$Q5.2 <= 5 | is.na(data$Q5.2)) & (data$Q5.2A <= 5 | is.na(data$Q5.2A)) & (data$Q5.2B <= 5 | is.na(data$Q5.2B)) &
             (data$Q5.3 <= 5 | is.na(data$Q5.3)) & (data$Q5.3A <= 5 | is.na(data$Q5.3A)) & (data$Q5.3B <= 5 | is.na(data$Q5.3B)) &
             (data$Q5.4 <= 5 | is.na(data$Q5.4)) & (data$Q5.4A <= 5 | is.na(data$Q5.4A)) & (data$Q5.4B <= 5 | is.na(data$Q5.4B)), ]

data.q1 <- data[(data$Q5.1 <= 5 | is.na(data$Q5.1)) & (data$Q5.1A <= 5 | is.na(data$Q5.1A)) & (data$Q5.1B <= 5 | is.na(data$Q5.1B)) , ]

data.q2 <- data[(data$Q5.2 <= 5 | is.na(data$Q5.2)) & (data$Q5.2A <= 5 | is.na(data$Q5.2A)) & (data$Q5.2B <= 5 | is.na(data$Q5.2B)) , ]

data.q3 <- data[ (data$Q5.3 <= 5 | is.na(data$Q5.3)) & (data$Q5.3A <= 5 | is.na(data$Q5.3A)) & (data$Q5.3B <= 5 | is.na(data$Q5.3B))  , ]

data.q4 <- data[ (data$Q5.4 <= 5 | is.na(data$Q5.4)) & (data$Q5.4A <= 5 | is.na(data$Q5.4A)) & (data$Q5.4B <= 5 | is.na(data$Q5.4B)) , ]

list.ISAF <- (data.mean$Q4C[data.mean$treat=="ISAF"] )

list.ISAF.q1 <- (data.q1$Q4C[data.q1$treat=="ISAF"] )

list.ISAF.q2 <- (data.q2$Q4C[data.q2$treat=="ISAF"] )

list.ISAF.q3 <- (data.q3$Q4C[data.q3$treat=="ISAF"] )

list.ISAF.q4 <- (data.q4$Q4C[data.q4$treat=="ISAF"] )

endorse1.ISAF <- (data.q1$Q5.1A[data.q1$treat=="ISAF"]) 

endorse2.ISAF <- (data.q2$Q5.2A[data.q2$treat=="ISAF"])

endorse3.ISAF <- (data.q3$Q5.3A[data.q3$treat=="ISAF"])

endorse4.ISAF <- (data.q4$Q5.4A[data.q4$treat=="ISAF"])

mean.ISAF <- apply(cbind(data.mean$Q5.1A[data.mean$treat=="ISAF"], data.mean$Q5.2A[data.mean$treat=="ISAF"],
                         data.mean$Q5.3A[data.mean$treat=="ISAF"], data.mean$Q5.4A[data.mean$treat=="ISAF"]), 1, mean)

list.control <- (data.mean$Q4A[data.mean$treat=="control"] )

list.control.q1 <- (data.q1$Q4A[data.q1$treat=="control"] )

list.control.q2 <- (data.q2$Q4A[data.q2$treat=="control"] )

list.control.q3 <- (data.q3$Q4A[data.q3$treat=="control"] )

list.control.q4 <- (data.q4$Q4A[data.q4$treat=="control"] )

endorse1.control <- (data.q1$Q5.1[data.q1$treat=="control"]) 

endorse2.control <- (data.q2$Q5.2[data.q2$treat=="control"])

endorse3.control <- (data.q3$Q5.3[data.q3$treat=="control"])

endorse4.control <- (data.q4$Q5.4[data.q4$treat=="control"])

mean.control <- apply(cbind(data.mean$Q5.1[data.mean$treat=="control"], data.mean$Q5.2[data.mean$treat=="control"],
                         data.mean$Q5.3[data.mean$treat=="control"], data.mean$Q5.4[data.mean$treat=="control"]), 1, mean)


###
###  Fig 2 - Comparison of Responses from List and Endorsement
##           Experiments in the Afghanistan Survey
###

pdf(file = "figs/CompareScatterEndorseListMean-ISAF.pdf", width = 8, height = 4)

par(mfcol = c(1,2), mar = c(4, 4, 1, 0), oma = c(0, 0, 0, 0))
u.control <- as.data.frame(unique(cbind(mean.control, list.control)))
u.control$n <- NA
for(i in 1:nrow(u.control))
  u.control$n[i] <- sum(mean.control == u.control[i,1] & list.control == u.control[i,2])

plot(u.control$mean.control, u.control$list.control, cex = u.control$n / length(mean.control)*60,
     ylim = c(-1, 5), xlim = c(0,6), main = "Control Group", xlab = "Endorsement Experiment",
     ylab = "List Experiment", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(mean.control, predict(loess(list.control ~ mean.control)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(mean.control, list.control),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(mean.control, list.control, method = "kendall"),2), nsmall = 2)), cex = 0.8)

u.treat <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u.treat$n <- NA
for(i in 1:nrow(u.treat))
  u.treat$n[i] <- sum(mean.ISAF == u.treat[i,1] & list.ISAF == u.treat[i,2])

plot(u.treat$mean.ISAF, u.treat$list.ISAF, cex = u.treat$n / length(mean.ISAF)*60,
     ylim = c(-1, 5), xlim = c(0,6), main = "ISAF Treatment Group",
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(mean.ISAF, list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(mean.ISAF, list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

dev.off()
comp.p <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = FALSE)
print(paste(c("overall", round(comp.p$ci, 3)), collapse = " "))

comp.p.kendall <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = FALSE, method = "kendall")
print(paste(c("overall tau", round(comp.p.kendall$ci, 3)), collapse = " "))


###
###  Fig 3 - Comparison of Responses from List and Endorsement
###          Experiments in the Afghanistan Survey by Endorsement Question
###

pdf(file = "figs/CompareScatterEndorseList-ISAF.pdf", width = 12, height = 6)
par(mfcol = c(2,4), mar = c(4, 4, 0, 0), oma = c(0, 1.15, 2, 0), cex = 1)
  
u <- as.data.frame(unique(cbind(endorse1.ISAF, list.ISAF.q1)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse1.ISAF == u[i,1] & list.ISAF.q1 == u[i,2])

plot(u$endorse1.ISAF, u$list.ISAF.q1, cex = u$n / length(endorse1.ISAF)*60, ylim = c(-1, 5), xlim = c(0,6), xlab = "", ylab = "List Experiment", axes = F, cex.main = 0.9)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse1.ISAF, predict(loess(list.ISAF.q1 ~ endorse1.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse1.ISAF, list.ISAF.q1),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse1.ISAF, list.ISAF.q1, method = "kendall"),2), nsmall = 2)), cex = 0.8)

comp.p <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",1, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE)
print(paste(c("direct elections", round(comp.p$ci, 3)), collapse = " "))

comp.p.kendall <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",1, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE, method = "kendall")
print(paste(c("direct elections tau", round(comp.p.kendall$ci, 3)), collapse = " "))

mtext(side = 3, paste("Direct Elections (p < .01)", sep = ""), line = 0.5 - .8, cex = .8, font = 2)

mtext(side = 2, "ISAF Treatment Group", line = 4.3, font = 2)
mtext(side = 2, "Control Group", line = 4.3, at = -7, font = 2)

u <- as.data.frame(unique(cbind(endorse1.control, list.control.q1)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse1.control == u[i,1] & list.control.q1 == u[i,2])

plot(u$endorse1.control, u$list.control.q1, cex = u$n / length(endorse1.control)*60, ylim = c(-1, 5), xlim = c(0,6),
     xlab = "Endorsement Experiment", ylab = "List Experiment", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse1.control, predict(loess(list.control.q1 ~ endorse1.control)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse1.control, list.control.q1),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse1.control, list.control.q1, method = "kendall"),2), nsmall = 2)), cex = 0.8)

u <- as.data.frame(unique(cbind(endorse2.ISAF, list.ISAF.q2)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse2.ISAF == u[i,1] & list.ISAF.q2 == u[i,2])

plot(u$endorse2.ISAF, u$list.ISAF.q2, cex = u$n / length(endorse2.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),  xlab = "", ylab = "", axes = F, cex.main = 0.9)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse2.ISAF, predict(loess(list.ISAF.q2 ~ endorse2.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse2.ISAF, list.ISAF.q2),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse2.ISAF, list.ISAF.q2, method = "kendall"),2), nsmall = 2)), cex = 0.8)

u <- as.data.frame(unique(cbind(endorse2.control, list.control.q2)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse2.control == u[i,1] & list.control.q2 == u[i,2])

comp.p <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",2, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE)
print(paste(c("prison reform", round(comp.p$ci, 3)), collapse = " "))

comp.p.kendall <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",2, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE, method = "kendall")
print(paste(c("prison reform tau", round(comp.p.kendall$ci, 3)), collapse = " "))

mtext(side = 3, paste("Prison Reform (p = ", format(round(comp.p$p.value, 2), nsmall = 2), ")", sep = ""), line = 0.5 - .8, cex = .8, font = 2)

plot(u$endorse2.control, u$list.control.q2, cex = u$n / length(endorse2.control)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse2.control, predict(loess(list.control.q2 ~ endorse2.control)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse2.control, list.control.q2),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse2.control, list.control.q2, method = "kendall"),2), nsmall = 2)), cex = 0.8)

u <- as.data.frame(unique(cbind(endorse3.ISAF, list.ISAF.q3)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse3.ISAF == u[i,1] & list.ISAF.q3 == u[i,2])

plot(u$endorse3.ISAF, u$list.ISAF.q3, cex = u$n / length(endorse3.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),   xlab = "", ylab = "", axes = F, cex.main = 0.9)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse3.ISAF, predict(loess(list.ISAF.q3 ~ endorse3.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse3.ISAF, list.ISAF.q3),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse3.ISAF, list.ISAF.q3, method = "kendall"),2), nsmall = 2)), cex = 0.8)

comp.p <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",3, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE)
print(paste(c("elec comm", round(comp.p$ci, 3)), collapse = " "))
mtext(side = 3, paste("Election Commission (p < .01)", sep = ""), line = 0.5 - .8, cex = .8, font = 2)

comp.p.kendall <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",3, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE, method = "kendall")
print(paste(c("elec comm tau", round(comp.p.kendall$ci, 3)), collapse = " "))


u <- as.data.frame(unique(cbind(endorse3.control, list.control.q3)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse3.control == u[i,1] & list.control.q3 == u[i,2])

plot(u$endorse3.control, u$list.control.q3, cex = u$n / length(endorse3.control)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),  ##main = paste("Election Commission (N = ", sum(u$n), ")", sep = ""),
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse3.control, predict(loess(list.control.q3 ~ endorse3.control)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse3.control, list.control.q3),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse3.control, list.control.q3, method = "kendall"),2), nsmall = 2)), cex = 0.8)

u <- as.data.frame(unique(cbind(endorse4.ISAF, list.ISAF.q4)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse4.ISAF == u[i,1] & list.ISAF.q4 == u[i,2])

plot(u$endorse4.ISAF, u$list.ISAF.q4, cex = u$n / length(endorse4.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),  xlab = "", ylab = "", axes = F, cex.main = 0.9)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse4.ISAF, predict(loess(list.ISAF.q4 ~ endorse4.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

comp.p <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",4, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE)
print(paste(c("corruption", round(comp.p$ci, 3)), collapse = " "))

comp.p.kendall <- comp.listEndorse(y.endorse = cor.dat[,paste("Q",4, sep = "")],
                    y.list = cor.dat[,"list.y"], treat = cor.dat[,"treat"],
                              endorse.mean = TRUE, method = "kendall")
print(paste(c("corruption tau", round(comp.p.kendall$ci, 3)), collapse = " "))
            
mtext(side = 3, paste("Corruption Reform (p < .01)", sep = ""), line = 0.5 - .8, cex = .8, font = 2)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse4.ISAF, list.ISAF.q4),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse4.ISAF, list.ISAF.q4,method = "kendall"),2), nsmall = 2)), cex = 0.8)

u <- as.data.frame(unique(cbind(endorse4.control, list.control.q4)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(endorse4.control == u[i,1] & list.control.q4 == u[i,2])

plot(u$endorse4.control, u$list.control.q4, cex = u$n / length(endorse4.control)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

pred <- cbind(endorse4.control, predict(loess(list.control.q4 ~ endorse4.control)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red") 

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(endorse4.control, list.control.q4),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(endorse4.control, list.control.q4, method = "kendall"),2), nsmall = 2)), cex = 0.8)

dev.off()

##
## Now create the subset figures
##

data <- data[(data$Q5.1 <= 5 | is.na(data$Q5.1)) & (data$Q5.1A <= 5 | is.na(data$Q5.1A)) & (data$Q5.1B <= 5 | is.na(data$Q5.1B)) &
             (data$Q5.2 <= 5 | is.na(data$Q5.2)) & (data$Q5.2A <= 5 | is.na(data$Q5.2A)) & (data$Q5.2B <= 5 | is.na(data$Q5.2B)) &
             (data$Q5.3 <= 5 | is.na(data$Q5.3)) & (data$Q5.3A <= 5 | is.na(data$Q5.3A)) & (data$Q5.3B <= 5 | is.na(data$Q5.3B)) &
             (data$Q5.4 <= 5 | is.na(data$Q5.4)) & (data$Q5.4A <= 5 | is.na(data$Q5.4A)) & (data$Q5.4B <= 5 | is.na(data$Q5.4B)), ]

data$violent.prov <- data$district.violence.count.5km >= mean(district.violence.count.5km$district.violence.count.5km)

###
###  Fig 4 - Comparison of Responses from a List Experiment and
###  Four Endorsement Experiments for ISAF by Levels of Violence and Territorial Control.
###

pdf(file = paste("figs/CompareScatterByViolenceControl-ISAF.pdf", sep=""), width = 12, height = 6)

par(mfrow = c(2,4), mar = c(4, 4, 0, 0), oma = c(0, 1.15, 2, 0), cex = 1)
  
## violent

## set up data

dat <- data[data$violent.prov == FALSE, ]

comp.p <- comp.listEndorse(y.endorse = cor.dat[data$violent.prov == FALSE,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[data$violent.prov == FALSE,"list.y"],
                              treat = cor.dat[data$violent.prov == FALSE,"treat"],
                              endorse.mean = FALSE)

print(paste(c("violence low", round(comp.p$ci, 3)), collapse = " "))

dataset <- 1

list.ISAF <- (dat$Q4C[dat$treat=="ISAF"] )
     
endorse1.ISAF <- (dat$Q5.1A[dat$treat=="ISAF"] )

endorse2.ISAF <- (dat$Q5.2A[dat$treat=="ISAF"])

endorse3.ISAF <- (dat$Q5.3A[dat$treat=="ISAF"] )

endorse4.ISAF <- (dat$Q5.4A[dat$treat=="ISAF"])

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "", ylab = "List Experiment", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

level.str <- c("District","Village")

mtext(side = 3, paste("Low (p < .01)", sep = ""), line = 0.5 - .8, cex = .8)

mtext(side = 3, paste(as.character(level.str[dataset]), "Violence Level"), line = 2 - .8, cex = 0.9, at = 7.3, font = 2)

mtext(side = 2, "ISAF Treatment Group", line = 4.3, font = 2)
mtext(side = 2, "Control Group", line = 4.3, at = -6.3, font = 2)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## non violent

## set up data

dat <- data[data$violent.prov == TRUE, ]

comp.p <- comp.listEndorse(y.endorse = cor.dat[data$violent.prov == TRUE,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[data$violent.prov == TRUE,"list.y"],
                              treat = cor.dat[data$violent.prov == TRUE,"treat"],
                              endorse.mean = FALSE)

print(paste(c("violence high", round(comp.p$ci, 3)), collapse = " "))

list.ISAF <- (dat$Q4C[dat$treat=="ISAF"] )
           
endorse1.ISAF <- (dat$Q5.1A[dat$treat=="ISAF"])

endorse2.ISAF <- (dat$Q5.2A[dat$treat=="ISAF"] )

endorse3.ISAF <- (dat$Q5.3A[dat$treat=="ISAF"])

endorse4.ISAF <- (dat$Q5.4A[dat$treat=="ISAF"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),
     xlab = "", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

mtext(side = 3, paste("High (p < .01)", sep = ""), line = 0.5 - .8, cex = .8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## control

## set up data

dat <- data[data$control == 0, ]
comp.p <- comp.listEndorse(y.endorse = cor.dat[data$control == 0,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[data$control == 0,"list.y"],
                              treat = cor.dat[data$control == 0,"treat"],
                              endorse.mean = FALSE)
print(paste(c("control taliban", round(comp.p$ci, 3)), collapse = " "))

list.ISAF <- (dat$Q4C[dat$treat=="ISAF"] )
     
endorse1.ISAF <- (dat$Q5.1A[dat$treat=="ISAF"] )

endorse2.ISAF <- (dat$Q5.2A[dat$treat=="ISAF"] )

endorse3.ISAF <- (dat$Q5.3A[dat$treat=="ISAF"] )

endorse4.ISAF <- (dat$Q5.4A[dat$treat=="ISAF"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

mtext(side = 3, paste("Taliban (p < .01)", sep = ""), line = 0.5 - .8, cex = .8)
mtext(side = 3, paste("District Territory Control"), line = 2 - .8, cex = 0.9, at = 7.3, font = 2)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## non violent

## set up data

dat <- data[data$control == 1, ]
comp.p <- comp.listEndorse(y.endorse = cor.dat[data$control == 1,paste("Q",1:4, sep = "")],
                    y.list = cor.dat[data$control == 1,"list.y"],
                              treat = cor.dat[data$control == 1,"treat"],
                              endorse.mean = FALSE)
print(paste(c("control contested", round(comp.p$ci, 3)), collapse = " "))

list.ISAF <- (dat$Q4C[dat$treat=="ISAF"])
           
endorse1.ISAF <- (dat$Q5.1A[dat$treat=="ISAF"])

endorse2.ISAF <- (dat$Q5.2A[dat$treat=="ISAF"] )

endorse3.ISAF <- (dat$Q5.3A[dat$treat=="ISAF"] )

endorse4.ISAF <- (dat$Q5.4A[dat$treat=="ISAF"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

mtext(side = 3, paste("Contested (p < .01)", sep = ""), line = -.3, cex = .8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## control group
dat <- data[data$violent.prov == FALSE, ]

list.ISAF <- (dat$Q4A[dat$treat=="control"] )
     
endorse1.ISAF <- (dat$Q5.1[dat$treat=="control"] )

endorse2.ISAF <- (dat$Q5.2[dat$treat=="control"])

endorse3.ISAF <- (dat$Q5.3[dat$treat=="control"] )

endorse4.ISAF <- (dat$Q5.4[dat$treat=="control"])

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "Endorsement Experiment", ylab = "List Experiment", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## non violent

## set up data

dat <- data[data$violent.prov == TRUE, ]

list.ISAF <- (dat$Q4A[dat$treat=="control"] )
           
endorse1.ISAF <- (dat$Q5.1[dat$treat=="control"])

endorse2.ISAF <- (dat$Q5.2[dat$treat=="control"] )

endorse3.ISAF <- (dat$Q5.3[dat$treat=="control"])

endorse4.ISAF <- (dat$Q5.4[dat$treat=="control"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## control

## set up data

dat <- data[data$control == 0, ]

list.ISAF <- (dat$Q4A[dat$treat=="control"] )
     
endorse1.ISAF <- (dat$Q5.1[dat$treat=="control"] )

endorse2.ISAF <- (dat$Q5.2[dat$treat=="control"] )

endorse3.ISAF <- (dat$Q5.3[dat$treat=="control"] )

endorse4.ISAF <- (dat$Q5.4[dat$treat=="control"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75), 
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")

## non violent

## set up data

dat <- data[data$control == 1, ]

list.ISAF <- (dat$Q4A[dat$treat=="control"])
           
endorse1.ISAF <- (dat$Q5.1[dat$treat=="control"])

endorse2.ISAF <- (dat$Q5.2[dat$treat=="control"] )

endorse3.ISAF <- (dat$Q5.3[dat$treat=="control"] )

endorse4.ISAF <- (dat$Q5.4[dat$treat=="control"] )

mean.ISAF <- apply(cbind(endorse1.ISAF, endorse2.ISAF, endorse3.ISAF, endorse4.ISAF), 1, mean)

## plot

u <- as.data.frame(unique(cbind(mean.ISAF, list.ISAF)))
u$n <- NA
for(i in 1:nrow(u))
  u$n[i] <- sum(mean.ISAF == u[i,1] & list.ISAF == u[i,2])

plot(u$mean.ISAF, u$list.ISAF, cex = u$n / length(mean.ISAF)*60,  ylim = c(-.75, 4.75), xlim = c(0.25,5.75),
     xlab = "Endorsement Experiment", ylab = "", axes = F)

axis(2, at = 0:4)
axis(1, at = 1:5)

text(1.5, 4.25, expression(rho), cex = 0.8)
text(1.5, 4.25, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF),2), nsmall = 2)), cex = 0.8)
text(1.5, 3.75, expression(tau), cex = 0.8)
text(1.5, 3.75, pos = 4, paste("= ", format(round(cor(u$mean.ISAF, u$list.ISAF, method = "kendall"),2), nsmall = 2)), cex = 0.8)

pred <- cbind(mean.ISAF, predict(loess(list.ISAF ~ mean.ISAF)))
pred <- pred[order(pred[,1]),]
lines(unique(pred), col = "red")


dev.off()

##
## Now do the regressions

library(arm)

load("results-list.RData")

res.delta.list.villdist <- list()
for(i in 1:length(res))
  res.delta.list.villdist[[i]] <- as.mcmc(res[[i]]$delta[(nrow(res[[i]]$delta)/2):nrow(res[[i]]$delta),])
res.delta.list.villdist <- as.mcmc.list(res.delta.list.villdist)

delta.est.list.villdist <- do.call(rbind, res.delta.list.villdist)

load("results-endorse.RData")

res.lambda.endorse.villdist <- list()
for(i in 1:length(res))
  res.lambda.endorse.villdist[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)/2):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)/2):nrow(res[[i]]$omega2),]))
res.lambda.endorse.villdist <- as.mcmc.list(res.lambda.endorse.villdist)

lambda.est.endorse.villdist <- do.call(rbind, res.lambda.endorse.villdist)

load("results-combine.RData")

res.lambda.combine.villdist <- list()
for(i in 1:length(res))
  res.lambda.combine.villdist[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)/2):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)/2):nrow(res[[i]]$omega2),]))
res.lambda.combine.villdist <- as.mcmc.list(res.lambda.combine.villdist)

lambda.est.combine.villdist <- do.call(rbind, res.lambda.combine.villdist)

logistic <- function(x) exp(x)/(1+exp(x))

coef.names <- c("harm.taliban","harm.taliban.NA",
                        "harm.ISAF","harm.ISAF.NA","approach.taliban",
                        "approach.taliban.NA","approach.ISAF","approach.ISAF.NA",
                        "frequency.ISAF","educ.years","age.tens","income","income.NA",
                        "madrassa","pro.taliban","pro.taliban.NA") ##

coef.names.villdist <- c("harm.taliban","harm.taliban.NA",
                        "harm.ISAF","harm.ISAF.NA","approach.taliban",
                        "approach.taliban.NA","approach.ISAF","approach.ISAF.NA",
                        "frequency.ISAF","educ.years","age.tens","income","income.NA",
                        "madrassa","pro.taliban","pro.taliban.NA", colnames(villcov), colnames(distcov))

coefs.list.villdist <- round(cbind(apply(delta.est.list.villdist, 2, mean),
                                   apply(delta.est.list.villdist, 2, sd),
                       apply(delta.est.list.villdist, 2, function(x) quantile(x, .025)),
                       apply(delta.est.list.villdist, 2, function(x) quantile(x, .975))), 3)
rownames(coefs.list.villdist) <- c("Intercept",coef.names.villdist)
colnames(coefs.list.villdist) <- c("est","se","lower.ci","upper.ci")

coefs.endorse.villdist <- round(cbind(apply(lambda.est.endorse.villdist, 2, mean)[1:29],
                                      apply(lambda.est.endorse.villdist, 2, sd)[1:29],
                       apply(lambda.est.endorse.villdist, 2, function(x) quantile(x, .025))[1:29],
                       apply(lambda.est.endorse.villdist, 2, function(x) quantile(x, .975))[1:29]), 3)
rownames(coefs.endorse.villdist) <- c("Intercept",coef.names.villdist)
colnames(coefs.endorse.villdist) <- c("est","se","lower.ci","upper.ci")

coefs.combine.villdist <- round(cbind(apply(lambda.est.combine.villdist, 2, mean),
                                      apply(lambda.est.combine.villdist, 2, sd),
                                      apply(lambda.est.combine.villdist, 2, function(x) quantile(x, .025)),
                                      apply(lambda.est.combine.villdist, 2, function(x) quantile(x, .975))), 3)
rownames(coefs.combine.villdist) <- c("Intercept",coef.names.villdist)
colnames(coefs.combine.villdist) <- c("est","se","lower.ci","upper.ci")

predict.logistic <- function(x.dat, x.post) logistic( as.matrix(x.dat) %*%
                           t(x.post[seq(from = 1, to = nrow(x.post), by = 5), 1:ncol(x.dat)]))

predict.probit <- function(x.dat, x.post) pnorm( as.matrix(x.dat) %*%
                           t(x.post[seq(from = 1, to = nrow(x.post), by = 5), 1:ncol(x.dat)]))

summ <- function(x) c(mean(x), quantile(x, .025), quantile(x, .975))
summ.diff <- function(x, y) c(mean(x - y), quantile(x - y, .025), quantile(x - y, .975))

predict.all <- function(x.data, x.post, func = predict.logistic) {

  harm.taliban.dat <- noharm.taliban.dat <- harm.ISAF.dat <- noharm.ISAF.dat <-
    approach.taliban.dat <- noapproach.taliban.dat <- approach.ISAF.dat <- noapproach.ISAF.dat <- overall.dat <- 
      as.data.frame(x.data)
  
  harm.taliban.dat$violent.exp.taliban <- 1
  noharm.taliban.dat$violent.exp.taliban <- 0
  harm.ISAF.dat$violent.exp.ISAF <- 1
  noharm.ISAF.dat$violent.exp.ISAF <- 0
  
  approach.taliban.dat$encounter.taliban <- 1
  noapproach.taliban.dat$encounter.taliban <- 0
  approach.ISAF.dat$encounter.ISAF <- 1
  noapproach.ISAF.dat$encounter.ISAF <- 0
  
  harm.taliban.pred <- apply(func(harm.taliban.dat, x.post), 2, mean)
  noharm.taliban.pred <- apply(func(noharm.taliban.dat, x.post), 2, mean)
  harm.ISAF.pred <- apply(func(harm.ISAF.dat, x.post), 2, mean)
  noharm.ISAF.pred <- apply(func(noharm.ISAF.dat, x.post), 2, mean)
  approach.taliban.pred <- apply(func(approach.taliban.dat, x.post), 2, mean)
  noapproach.taliban.pred <- apply(func(noapproach.taliban.dat, x.post), 2, mean)
  approach.ISAF.pred <- apply(func(approach.ISAF.dat, x.post), 2, mean)
  noapproach.ISAF.pred <- apply(func(noapproach.ISAF.dat, x.post), 2, mean)
  mean.pred <- apply(func(overall.dat, x.post), 2, mean)

  return.mat <- rbind(summ(harm.taliban.pred), summ(noharm.taliban.pred),
                      summ.diff(harm.taliban.pred, noharm.taliban.pred),
                      summ(harm.ISAF.pred), summ(noharm.ISAF.pred),
                      summ.diff(harm.ISAF.pred, noharm.ISAF.pred),
                      summ(approach.taliban.pred), summ(noapproach.taliban.pred),
                      summ.diff(approach.taliban.pred, noapproach.taliban.pred),
                      summ(approach.ISAF.pred), summ(noapproach.ISAF.pred),
                      summ.diff(approach.ISAF.pred, noapproach.ISAF.pred),
                      summ(mean.pred))

  rownames(return.mat) <- c("harm taliban","no harm taliban", "diff taliban (harm - no harm)",
                            "harm ISAF","no harm ISAF", "diff ISAF (harm - no harm)",
                            "approach taliban","no approach taliban", "diff taliban (approach - no approach)",
                            "approach ISAF","no approach ISAF", "diff ISAF (approach - no approach)", "mean overall")

  return(return.mat)
  
}

predict.diff <- function(x.data, x.post.list, x.post.endorse,
                        func = predict.logistic) {

  harm.taliban.dat <- noharm.taliban.dat <- harm.ISAF.dat <- noharm.ISAF.dat <-
    approach.taliban.dat <- noapproach.taliban.dat <- approach.ISAF.dat <- noapproach.ISAF.dat <- overall.dat <- 
      as.data.frame(x.data)
  
  harm.taliban.dat$violent.exp.taliban <- 1
  noharm.taliban.dat$violent.exp.taliban <- 0
  harm.ISAF.dat$violent.exp.ISAF <- 1
  noharm.ISAF.dat$violent.exp.ISAF <- 0
  
  approach.taliban.dat$encounter.taliban <- 1
  noapproach.taliban.dat$encounter.taliban <- 0
  approach.ISAF.dat$encounter.ISAF <- 1
  noapproach.ISAF.dat$encounter.ISAF <- 0

  n.draws <- min(nrow(x.post.list), nrow(x.post.endorse))
  
  harm.taliban.pred <- apply(func(harm.taliban.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(harm.taliban.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  noharm.taliban.pred <- apply(func(noharm.taliban.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(noharm.taliban.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  harm.ISAF.pred <- apply(func(harm.ISAF.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(harm.ISAF.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  noharm.ISAF.pred <- apply(func(noharm.ISAF.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(noharm.ISAF.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  approach.taliban.pred <- apply(func(approach.taliban.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(approach.taliban.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  noapproach.taliban.pred <- apply(func(noapproach.taliban.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(noapproach.taliban.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  approach.ISAF.pred <- apply(func(approach.ISAF.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(approach.ISAF.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  noapproach.ISAF.pred <- apply(func(noapproach.ISAF.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(noapproach.ISAF.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  mean.pred <- apply(func(overall.dat,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(overall.dat,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)

  return.mat <- rbind(summ(harm.taliban.pred), summ(noharm.taliban.pred),
                      summ.diff(harm.taliban.pred, noharm.taliban.pred),
                      summ(harm.ISAF.pred), summ(noharm.ISAF.pred),
                      summ.diff(harm.ISAF.pred, noharm.ISAF.pred),
                      summ(approach.taliban.pred), summ(noapproach.taliban.pred),
                      summ.diff(approach.taliban.pred, noapproach.taliban.pred),
                      summ(approach.ISAF.pred), summ(noapproach.ISAF.pred),
                      summ.diff(approach.ISAF.pred, noapproach.ISAF.pred),
                      summ(mean.pred))

  rownames(return.mat) <- c("harm taliban","no harm taliban", "diff taliban (harm - no harm)",
                            "harm ISAF","no harm ISAF", "diff ISAF (harm - no harm)",
                            "approach taliban","no approach taliban", "diff taliban (approach - no approach)",
                            "approach ISAF","no approach ISAF", "diff ISAF (approach - no approach)", "mean overall")

  return(return.mat)
  
}

list.villdist.pred <- predict.all(villdist.data, delta.est.list.villdist,
                                  predict.probit)

endorse.villdist.pred <- predict.all(villdist.data, lambda.est.endorse.villdist,
                                     predict.probit)

combine.villdist.pred <- predict.all(villdist.data, lambda.est.combine.villdist,
                                     predict.probit)

diff.villdist.pred <- predict.diff(villdist.data, delta.est.list.villdist,
                                  lambda.est.endorse.villdist,
                                  predict.probit)

list.indiv.pred <- predict.probit(villdist.data, delta.est.list.villdist)
list.indiv.pred.means <- apply(list.indiv.pred, 1, mean)
list.indiv.pred.sd <- apply(list.indiv.pred, 1, sd)
endorse.indiv.pred <- predict.probit(villdist.data, lambda.est.endorse.villdist)
endorse.indiv.pred.means <- apply(endorse.indiv.pred, 1, mean)
endorse.indiv.pred.sd <- apply(endorse.indiv.pred, 1, sd)
diff.indiv.pred.means <- list.indiv.pred.means - endorse.indiv.pred.means

predict.setx <- function(x.data, var, setx = NULL, x.post.list, x.post.endorse, x.post.combine, func = predict.probit) {

  n.draws <- min(nrow(x.post.list), nrow(x.post.endorse))

  if(var != "control") {
    x.data[,paste(var)] <- setx
  } else {
    x.data[,paste("control.gov")] <- setx == 2
    x.data[,paste("control.contested")] <- setx == 3
  }
  
  pred.list <- apply(func(x.data, x.post.list), 2, mean)
  pred.endorse <- apply(func(x.data, x.post.endorse), 2, mean)
  pred.diff <- apply(func(x.data,
                                  x.post.list[(nrow(x.post.list)-n.draws):
                                              nrow(x.post.list), ]) -
                             func(x.data,
                                  x.post.endorse[(nrow(x.post.endorse)-n.draws):
                                              nrow(x.post.endorse), ]), 2, mean)
  pred.combine <- apply(func(x.data, x.post.combine), 2, mean)

  return(c(summ(pred.list), summ(pred.endorse), summ(pred.diff), summ(pred.combine)))
  
}

data$cerp.money <- data$cerp.money.std * 722365.5 + 608540.2

lower.bound <- ( quantile(data$cerp.money, .1) - 608540.2 ) / 722365.5
upper.bound <- ( quantile(data$cerp.money, .9) - 608540.2 ) / 722365.5

cerp.vals <- seq(from = lower.bound, to = upper.bound,
                 length.out = 15)
cerp.pred <- matrix(NA, ncol = 3*4, nrow = length(cerp.vals))
for(i in 1:length(cerp.vals))
  cerp.pred[i,] <- predict.setx(villdist.data, "cerp.money.std", cerp.vals[i], delta.est.list.villdist,
                                lambda.est.endorse.villdist, lambda.est.combine.villdist, predict.probit)
colnames(cerp.pred) <- rep(c("est","lower.ci","upper.ci"), 4)
rownames(cerp.pred) <- cerp.vals

control.vals <- c(1,2,3)
control.pred <- matrix(NA, ncol = 3*4, nrow = length(control.vals))
for(i in 1:length(control.vals))
  control.pred[i,] <- predict.setx(villdist.data, "control", control.vals[i], delta.est.list.villdist,
                                lambda.est.endorse.villdist, lambda.est.combine.villdist, predict.probit)
colnames(control.pred) <- rep(c("est","lower.ci","upper.ci"), 4)
rownames(control.pred) <- control.vals

##
## Fig 5 - Estimated Overall Proportion of ISAF Supporters based on the List Experiment, En- dorsement Experiment, and Combined Models
##

pdf("figs/SupportCompare.pdf", height = 3, width = 4)
par(mfcol = c(1,1), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.25, 1.5, 1.75, 2)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1.15, 2.1), ylim = c(-.25, 0.5), axes = FALSE)

axis(2, at = seq(from = -.25, to = .5, by = 0.25), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Overall Proportion of ISAF Supporters", outer = TRUE, side = 2, cex = 0.65,
      at = 0.5, line = 1.5, padj = .1)

coefsPlot <- rbind(list.villdist.pred["mean overall",],
                   endorse.villdist.pred["mean overall",],
                   diff.villdist.pred["mean overall",],
                   combine.villdist.pred["mean overall",]
                   )
colnames(coefsPlot) <- c("est", "lower.ci","upper.ci")

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 3), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

text(x.coords[1] , coefsPlot[2,"upper.ci"] + .2, "List", cex = 1.15)
text(x.coords[2], coefsPlot[2,"upper.ci"] + .2, "Endorsement", cex = 1.15)
text(x.coords[3], coefsPlot[2,"upper.ci"] + .213, "Difference\n(List - Endorse)", cex = 1.15, adj = c(0.5, 1))
text(x.coords[4], coefsPlot[2,"upper.ci"] + .2, "Combined", cex = 1.15)

dev.off()

##
## Figure 6 Panel A - Estimated Average Marginal Effects of Taliban and ISAF Victimization on the Probability of Supporting ISAF
## 

pdf("figs/CoefsHarmCompare.pdf", height = 3, width = 4)
par(mfcol = c(1,1), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.33, 1.66, 1.99, 2.99, 3.33, 3.66)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1, 3.99), ylim = c(-.25, .5), axes = FALSE)

axis(2, at = seq(from = -.25, to = .5, by = 0.25), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Effects on Probability of Supporting ISAF", outer = TRUE, side = 2, cex = 0.65,
      at = 0.5, line = 1.5, padj = .1)

coefsPlot <- rbind(list.villdist.pred["diff taliban (harm - no harm)",],
                   endorse.villdist.pred["diff taliban (harm - no harm)",],
                   combine.villdist.pred["diff taliban (harm - no harm)",],

                   list.villdist.pred["diff ISAF (harm - no harm)",],
                   endorse.villdist.pred["diff ISAF (harm - no harm)",],
                   combine.villdist.pred["diff ISAF (harm - no harm)",]
                   )
colnames(coefsPlot) <- c("est", "lower.ci","upper.ci")

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 3), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

text(x.coords[1] , coefsPlot[1,"lower.ci"] - .1, "List", cex = 1)
text(x.coords[2], coefsPlot[2,"upper.ci"] + .1, "Endorse", cex = 1)
text(x.coords[3], coefsPlot[3,"lower.ci"] - .1, "Combined", cex = 1)

text(x.coords[4] , coefsPlot[6,"lower.ci"] - .1, "List", cex = 1)
text(x.coords[5], coefsPlot[5,"upper.ci"] +  .1, "Endorse", cex = 1)
text(x.coords[6], coefsPlot[6,"lower.ci"] - .1, "Combined", cex = 1)

mtext("Victimization", outer = TRUE, cex = 0.7, at = 0.475, line = 0, padj = -1.65, font = 2)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.245, line = -.7)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.71, line = -.7)

dev.off()

##
## Figure 6 Panel B - Estimated Average Marginal Effects of Taliban and ISAF Post-harm Mitigation Efforts (“Approach” by Combatants) on the Probability of Supporting ISAF
## 

pdf("figs/CoefsApproachCompare.pdf", height = 3, width = 4)
par(mfcol = c(1,1), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.33, 1.66, 1.99, 2.99, 3.33, 3.66)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1, 3.99), ylim = c(-.25, .5), axes = FALSE)

axis(2, at = seq(from = -.25, to = .5, by = 0.25), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Effects on Probability of Supporting ISAF", outer = TRUE, side = 2, cex = 0.65,
      at = 0.5, line = 1.5, padj = .1)

coefsPlot <- rbind(list.villdist.pred["diff taliban (approach - no approach)",],
                   endorse.villdist.pred["diff taliban (approach - no approach)",],
                   combine.villdist.pred["diff taliban (approach - no approach)",],
                   
                   list.villdist.pred["diff ISAF (approach - no approach)",],
                   endorse.villdist.pred["diff ISAF (approach - no approach)",],
                   combine.villdist.pred["diff ISAF (approach - no approach)",]
                   )
colnames(coefsPlot) <- c("est", "lower.ci","upper.ci")

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 3), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

text(x.coords[1] , coefsPlot[1,"upper.ci"] + .1, "List", cex = 1)
text(x.coords[2], coefsPlot[2,"upper.ci"] + .2, "Endorse", cex = 1)
text(x.coords[3], coefsPlot[3,"upper.ci"] + .15, "Combined", cex = 1)

text(x.coords[4] , coefsPlot[4,"lower.ci"] - .1, "List", cex = 1)
text(x.coords[5], coefsPlot[5,"lower.ci"] -  .1, "Endorse", cex = 1)
text(x.coords[6], coefsPlot[6,"lower.ci"] - .1, "Combined", cex = 1)

mtext("Approach after Victimization", outer = TRUE, cex = 0.7, at = 0.475, line = 0, padj = -1.65, font = 2)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.245, line = -.7)
mtext("by ISAF", outer = TRUE, cex = 0.65, at = 0.71, line = -.7)

dev.off()

##
## Fig 7 Panel A - Estimated Proportion of ISAF Supporters based on the List Experiment, Endorsement Experiment, and Combined Models As a Function of by the Amount of Aid
##

pdf("figs/CoefsCERPCompare.pdf", height = 3, width = 5.5)
par(mfcol = c(1,1), las = 1, mar=c(2,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

list.est <- cerp.pred[,1:3]
endorse.est <- cerp.pred[,4:6]
diff.est <- cerp.pred[,7:9]
combine.est <- cerp.pred[,10:12]

x.coords <- as.numeric(rownames(list.est))

## CERP was standardized in data.R using mean = 608540.2 and sd = 722365.5

x.coords <- (x.coords * 722365.5 + 608540.2) / 100000

### ISAF Support

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(min(x.coords), max(x.coords)), ylim = c(0,0.8), axes = FALSE,
     xlab = "CERP Aid Spending (hundred thousands)")

axis(2, at = seq(from = 0, to = 0.5, by = 0.25), las = 1)
axis(1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Proportion of ISAF Supporters", outer = TRUE, side = 2, cex = 0.65,
      at = 0.4, line = 1.5, padj = .1)
mtext("CERP Aid Spending (hundred thousands)", outer = TRUE, side = 1, cex = 0.65,
      at = 0.5, line = 0.95, padj = .1)

points(x.coords, list.est[,"est"], pch = 19, cex = 1.5)
segments(x.coords, list.est[,c("lower.ci")], x.coords, list.est[,c("upper.ci")])

points(x.coords + .25, endorse.est[,"est"], pch = 15, cex = 1.5)
segments(x.coords + .25, endorse.est[,c("lower.ci")], x.coords + .25, endorse.est[,c("upper.ci")])

text(0.5, 0.35, "Endorsement\nExperiment", adj = c(0, NA))
text(8.5, 0.5, "List\nExperiment", adj = c(1, NA))

dev.off()

##
## Fig 7 Panel B - Estimated Proportion of ISAF Supporters based on the List Experiment, Endorsement Experiment, and Combined Models As a Function of by the Territorial Control
##


pdf("figs/CoefsControlCompare.pdf", height = 3, width = 5.5)
par(mfcol = c(1,1), las = 1, mar=c(2,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

list.est <- control.pred[,1:3]
endorse.est <- control.pred[,4:6]
diff.est <- control.pred[,7:9]
combine.est <- control.pred[,10:12]

x.coords <- as.numeric(rownames(list.est)) / 1.1

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(min(x.coords)-.3, max(x.coords)+.7), ylim = c(0,0.8), axes = FALSE)

print(min(diff.est[,"lower.ci"]))

axis(2, at = seq(from = 0, to = 0.5, by = 0.25), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Proportion of ISAF Supporters", outer = TRUE, side = 2, cex = 0.65,
      at = 0.4, line = 1.5, padj = .1)

points(x.coords, list.est[,"est"], pch = 19, cex = 1.5)
segments(x.coords, list.est[,c("lower.ci")], x.coords, list.est[,c("upper.ci")])

points(x.coords + .2, endorse.est[,"est"], pch = 15, cex = 1.5)
segments(x.coords + .2, endorse.est[,c("lower.ci")], x.coords + .2, endorse.est[,c("upper.ci")])

points(x.coords + .4, combine.est[,"est"], pch = 16, cex = 1.5)
segments(x.coords + .4, combine.est[,c("lower.ci")], x.coords + .4, combine.est[,c("upper.ci")])

text(x.coords[1] , list.est[1 ,"upper.ci"] + .15, "List", cex = 1)
text(x.coords[1] + .2, endorse.est[1 ,"upper.ci"] + .15, "Endorse", cex = 1)
text(x.coords[1] + .4, combine.est[1 ,"upper.ci"] + .1, "Combine", cex = 1)

text(x.coords[2] , list.est[ 2,"upper.ci"] + .1, "List", cex = 1)
text(x.coords[2] + .2, endorse.est[ 2,"upper.ci"] + .08, "Endorse", cex = 1)
text(x.coords[2] + .4, combine.est[2 ,"upper.ci"] + .08, "Combine", cex = 1)

text(x.coords[3] , list.est[3 ,"upper.ci"] + .05, "List", cex = 1)
text(x.coords[3] + .2, endorse.est[3 ,"upper.ci"] + .1, "Endorse", cex = 1)
text(x.coords[3] + .4, combine.est[3 ,"upper.ci"] + .15, "Combine", cex = 1)

mtext("Territorial Control", outer = TRUE, cex = 0.7, at = 0.475, line = -2.3, padj = -1.7, font = 2)
mtext("by Taliban", outer = TRUE, cex = 0.65, at = 0.22, line = -2.5)
mtext("by the Government", outer = TRUE, cex = 0.65, at = 0.4755, line = -2.5)
mtext("Contested", outer = TRUE, cex = 0.65, at = 0.73, line = -2.5)

dev.off()

##
## Table A.4 - Estimated Coefficients for the Three Models
##

library(xtable)

tab <- cbind(coefs.list.villdist[,1:2], coefs.endorse.villdist[,1:2], coefs.combine.villdist[,1:2])

rownames(tab) <- c("Intercept", "Harm from Taliban violence", "Harm from Taliban violence is NA",
                   "Harm from ISAF violence", "Harm from ISAF violence is NA",
                   "Approach by Taliban after Harm", "Approach by Taliban after Harm is NA",
                   "Approach by ISAF after Harm", "Approach by ISAF after Harm is NA",
                   "ISAF encounter frequency", "Years of education",
                   "Age (tens)", "Income (Afghanis)", "Income is NA",
                   "Schooled in madrassa", "Pro-Taliban tribe",
                   "Pro-Taliban tribe is NA", "Altitude (km)",
                   "Population", "ISAF-initiated violent events (within 5km)",
                   "Taliban-initiated violent events (within 5km)",
                   "Sha'ria courts", "CERP project spending",
                   "Opium cultivation (ha.)", "CDC project count",
                   "Road length (km)", "Pakistan border",
                   "Government territorial control","Contested territorial control")
                   
sink("figs/compRegResults.txt")
print(xtable(tab))
sink()

## drop policy analysis

load("results-endorse-dropQ1.RData")

res.lambda.endorse.dropQ1 <- list()
for(i in 1:length(res))
  res.lambda.endorse.dropQ1[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)/2):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)/2):nrow(res[[i]]$omega2),]))
res.lambda.endorse.dropQ1 <- as.mcmc.list(res.lambda.endorse.dropQ1)

lambda.est.endorse.dropQ1 <- do.call(rbind, res.lambda.endorse.dropQ1)

coefs.endorse.dropQ1 <- round(cbind(apply(lambda.est.endorse.dropQ1, 2, mean)[1:29],
                                      apply(lambda.est.endorse.dropQ1, 2, sd)[1:29],
                       apply(lambda.est.endorse.dropQ1, 2, function(x) quantile(x, .025))[1:29],
                       apply(lambda.est.endorse.dropQ1, 2, function(x) quantile(x, .975))[1:29]), 3)
colnames(coefs.endorse.dropQ1) <- c("est","se","lower.ci","upper.ci")

load("results-endorse-dropQ2.RData")

res.lambda.endorse.dropQ2 <- list()
for(i in 1:length(res))
  res.lambda.endorse.dropQ2[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)/2):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)/2):nrow(res[[i]]$omega2),]))
res.lambda.endorse.dropQ2 <- as.mcmc.list(res.lambda.endorse.dropQ2)

lambda.est.endorse.dropQ2 <- do.call(rbind, res.lambda.endorse.dropQ2)

coefs.endorse.dropQ2 <- round(cbind(apply(lambda.est.endorse.dropQ2, 2, mean)[1:29],
                                      apply(lambda.est.endorse.dropQ2, 2, sd)[1:29],
                       apply(lambda.est.endorse.dropQ2, 2, function(x) quantile(x, .025))[1:29],
                       apply(lambda.est.endorse.dropQ2, 2, function(x) quantile(x, .975))[1:29]), 3)
colnames(coefs.endorse.dropQ2) <- c("est","se","lower.ci","upper.ci")

load("results-endorse-dropQ3.RData")

res.lambda.endorse.dropQ3 <- list()
for(i in 1:length(res))
  res.lambda.endorse.dropQ3[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)/2):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)/2):nrow(res[[i]]$omega2),]))
res.lambda.endorse.dropQ3 <- as.mcmc.list(res.lambda.endorse.dropQ3)

lambda.est.endorse.dropQ3 <- do.call(rbind, res.lambda.endorse.dropQ3)

coefs.endorse.dropQ3 <- round(cbind(apply(lambda.est.endorse.dropQ3, 2, mean)[1:29],
                                      apply(lambda.est.endorse.dropQ3, 2, sd)[1:29],
                       apply(lambda.est.endorse.dropQ3, 2, function(x) quantile(x, .025))[1:29],
                       apply(lambda.est.endorse.dropQ3, 2, function(x) quantile(x, .975))[1:29]), 3)
colnames(coefs.endorse.dropQ3) <- c("est","se","lower.ci","upper.ci")

load("results-endorse-dropQ4.RData")

res.lambda.endorse.dropQ4 <- list()
for(i in 1:length(res))
  res.lambda.endorse.dropQ4[[i]] <- mcmc(res[[i]]$lambda[(nrow(res[[i]]$lambda)-4000):nrow(res[[i]]$lambda),] /
                                           sqrt(res[[i]]$omega2[(nrow(res[[i]]$omega2)-4000):nrow(res[[i]]$omega2),]))
res.lambda.endorse.dropQ4 <- as.mcmc.list(res.lambda.endorse.dropQ4)

lambda.est.endorse.dropQ4 <- do.call(rbind, res.lambda.endorse.dropQ4)

coefs.endorse.dropQ4 <- round(cbind(apply(lambda.est.endorse.dropQ4, 2, mean)[1:29],
                                      apply(lambda.est.endorse.dropQ4, 2, sd)[1:29],
                       apply(lambda.est.endorse.dropQ4, 2, function(x) quantile(x, .025))[1:29],
                       apply(lambda.est.endorse.dropQ4, 2, function(x) quantile(x, .975))[1:29]), 3)
colnames(coefs.endorse.dropQ4) <- c("est","se","lower.ci","upper.ci")

tab <- cbind(coefs.endorse.villdist[,1:2], coefs.endorse.dropQ1[,1:2], coefs.endorse.dropQ2[,1:2], coefs.endorse.dropQ3[,1:2],
             coefs.endorse.dropQ4[,1:2])

rownames(tab) <- c("Intercept", "Harm from Taliban violence", "Harm from Taliban violence is NA",
                   "Harm from ISAF violence", "Harm from ISAF violence is NA",
                   "Approach by Taliban after Harm", "Approach by Taliban after Harm is NA",
                   "Approach by ISAF after Harm", "Approach by ISAF after Harm is NA",
                   "ISAF encounter frequency", "Years of education",
                   "Age (tens)", "Income (Afghanis)", "Income is NA",
                   "Schooled in madrassa", "Pro-Taliban tribe",
                   "Pro-Taliban tribe is NA", "Altitude (km)",
                   "Population", "ISAF-initiated violent events (within 5km)",
                   "Taliban-initiated violent events (within 5km)",
                   "Sha'ria courts", "CERP project spending",
                   "Opium cultivation (ha.)", "CDC project count",
                   "Road length (km)", "Pakistan border",
                   "Government territorial control","Contested territorial control")
                   
sink("figs/endorseDropPolicyResults.txt")
print(xtable(tab))
sink()

endorse.villdist.pred <- predict.all(villdist.data, lambda.est.endorse.villdist,
                                     predict.probit)

endorse.villdist.pred.drop1 <- predict.all(villdist.data, lambda.est.endorse.dropQ1,
                                     predict.probit)

endorse.villdist.pred.drop2 <- predict.all(villdist.data, lambda.est.endorse.dropQ2,
                                     predict.probit)

endorse.villdist.pred.drop3 <- predict.all(villdist.data, lambda.est.endorse.dropQ3,
                                     predict.probit)

endorse.villdist.pred.drop4 <- predict.all(villdist.data, lambda.est.endorse.dropQ4,
                                     predict.probit)

pdf("figs/SupportCompareDropPolicy.pdf", height = 2.5, width = 4)
par(mfcol = c(1,1), las = 1, mar=c(0,2,0,3.5), oma=c(2,3,2,0), cex = 0.6)

x.coords <- c(1.25, 1.5, 1.75, 2, 2.25)

plot(0:11, seq(from = -2, to = 1.5, length = 12), type = "n",
     xlim = c(1.15, 2.35), ylim = c(0, 0.5), axes = FALSE)

axis(2, at = seq(from = 0, to = .5, by = 0.25), las = 1)
abline(0, 0, lty = "dashed")
par(las = 0)
mtext("Overall Proportion of ISAF Supporters", outer = TRUE, side = 2, cex = 0.65,
      at = 0.5, line = 1.5, padj = .1)

coefsPlot <- rbind(endorse.villdist.pred["mean overall",],
                   endorse.villdist.pred.drop1["mean overall",],
                   endorse.villdist.pred.drop2["mean overall",],
                   endorse.villdist.pred.drop3["mean overall",],
                   endorse.villdist.pred.drop4["mean overall",]
                   )
colnames(coefsPlot) <- c("est", "lower.ci","upper.ci")

points(x.coords, coefsPlot[,"est"], pch = rep(c(19,15,17), 3), cex = 1.5)
segments(x.coords, coefsPlot[,c("lower.ci")], x.coords, coefsPlot[,c("upper.ci")])

text(x.coords[1] , coefsPlot[2,"upper.ci"] + .3, "All Four\nPolicies", cex = 1.15)
text(x.coords[2], coefsPlot[2,"upper.ci"] + .3, "Drop\nDirect\nElections", cex = 1.15)
text(x.coords[3], coefsPlot[2,"upper.ci"] + .3, "Drop\nPrison\nReform", cex = 1.15)
text(x.coords[4], coefsPlot[2,"upper.ci"] + .3, "Drop\nElection\nCommission", cex = 1.15)
text(x.coords[5] , coefsPlot[2,"upper.ci"] + .3, "Drop\nCorruption\nReform", cex = 1.15)

dev.off()
